/*
 *  Java Information Dynamics Toolkit (JIDT)
 *  Copyright (C) 2012, Joseph T. Lizier
 *  
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *  
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *  
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

package infodynamics.utils;

import infodynamics.utils.MathsUtils;

import junit.framework.TestCase;

/**
 * Test functionality of the utility functions in MathsUtils
 * 
 * @author Joseph Lizier.
 *
 */
public class MathsUtilsTest extends TestCase {

	private static double OCTAVE_RESOLUTION = 0.00001;

	/** 
	 * Confirm that our erf() function is correct to 6 dp
	 */
	public void testErf() {
		// Compare against the results supplied by octave
		double[] octaveResults = {0.000000, 0.112463, 0.222703, 0.328627,
				0.428392, 0.520500, 0.603856, 0.677801, 0.742101, 0.796908,
				0.842701, 0.880205, 0.910314, 0.934008, 0.952285, 0.966105,
				0.976348, 0.983790, 0.989091, 0.992790, 0.995322, 0.997021,
				0.998137, 0.998857, 0.999311, 0.999593, 0.999764, 0.999866,
				0.999925, 0.999959, 0.999978, 0.999988, 0.999994, 0.999997,
				0.999998, 0.999999, 1.000000};
		
		for (int n=0; n<octaveResults.length; n++) {
			assertEquals(octaveResults[n], MathsUtils.erf(n*0.1), 0.000001);
		}
	}

	/** 
	 * Confirm that our chicdf() function is correct to 6 decimal places
	 */
	public void testChiCdf() {
		// Compare against the results supplied by octave
		
		// First for 1 degree of freedom (odd, base case)
		double[] octaveResults1degFree = {0.000000, 0.248170, 0.345279, 0.416118,
				0.472911, 0.520500, 0.561422, 0.597216, 0.628907, 0.657218,
				0.682689, 0.705734, 0.726678, 0.745787, 0.763276, 0.779329,
				0.794097, 0.807712, 0.820288, 0.831922, 0.842701, 0.852701,
				0.861989, 0.870626, 0.878665, 0.886154, 0.893136, 0.899652,
				0.905736, 0.911420, 0.916735, 0.921708, 0.926362, 0.930720,
				0.934804, 0.938631, 0.942220, 0.945588, 0.948747, 0.951714,
				0.954500, 0.957117, 0.959576, 0.961888, 0.964061, 0.966105,
				0.968028, 0.969837, 0.971540, 0.973143, 0.974653, 0.976074,
				0.977413, 0.978675, 0.979863, 0.980984, 0.982040, 0.983035,
				0.983974, 0.984859, 0.985694, 0.986482, 0.987225, 0.987926,
				0.988588, 0.989213, 0.989802, 0.990359, 0.990884, 0.991380,
				0.991849, 0.992292, 0.992710, 0.993105, 0.993478, 0.993830,
				0.994163, 0.994478, 0.994775, 0.995057, 0.995322, 0.995573,
				0.995811, 0.996035, 0.996248, 0.996449, 0.996638, 0.996818,
				0.996988, 0.997148, 0.997300, 0.997444, 0.997580, 0.997708,
				0.997830, 0.997945, 0.998054, 0.998157, 0.998255, 0.998347,
				0.998435};
		for (int n=0; n<octaveResults1degFree.length; n++) {
			assertEquals(octaveResults1degFree[n],
					MathsUtils.chiSquareCdf(n*0.1, 1), 0.000001);
		}
		
		// Next for 2 degrees of freedom (even, base case)
		double[] octaveResults2degFree = {0.000000, 0.048771, 0.095163, 0.139292,
				0.181269, 0.221199, 0.259182, 0.295312, 0.329680, 0.362372,
				0.393469, 0.423050, 0.451188, 0.477954, 0.503415, 0.527633,
				0.550671, 0.572585, 0.593430, 0.613259, 0.632121, 0.650062,
				0.667129, 0.683363, 0.698806, 0.713495, 0.727468, 0.740760,
				0.753403, 0.765430, 0.776870, 0.787752, 0.798103, 0.807950,
				0.817316, 0.826226, 0.834701, 0.842763, 0.850431, 0.857726,
				0.864665, 0.871265, 0.877544, 0.883516, 0.889197, 0.894601,
				0.899741, 0.904631, 0.909282, 0.913706, 0.917915, 0.921918,
				0.925726, 0.929349, 0.932794, 0.936072, 0.939190, 0.942156,
				0.944977, 0.947660, 0.950213, 0.952641, 0.954951, 0.957148,
				0.959238, 0.961226, 0.963117, 0.964916, 0.966627, 0.968254,
				0.969803, 0.971275, 0.972676, 0.974009, 0.975276, 0.976482,
				0.977629, 0.978720, 0.979758, 0.980745, 0.981684, 0.982578,
				0.983427, 0.984236, 0.985004, 0.985736, 0.986431, 0.987093,
				0.987723, 0.988321, 0.988891, 0.989433, 0.989948, 0.990438,
				0.990905, 0.991348, 0.991770, 0.992172, 0.992553, 0.992917,
				0.993262};
		for (int n=0; n<octaveResults2degFree.length; n++) {
			assertEquals(octaveResults2degFree[n],
					MathsUtils.chiSquareCdf(n*0.1, 2), 0.000001);
		}

		// Next for 3 degrees of freedom (odd, > 1)
		double[] octaveResults3degFree = {0.000000, 0.008163, 0.022411, 0.039972,
				0.059758, 0.081109, 0.103568, 0.126796, 0.150533, 0.174572,
				0.198748, 0.222926, 0.246996, 0.270867, 0.294465, 0.317730,
				0.340610, 0.363066, 0.385065, 0.406581, 0.427593, 0.448087,
				0.468052, 0.487479, 0.506365, 0.524709, 0.542510, 0.559773,
				0.576500, 0.592698, 0.608375, 0.623537, 0.638195, 0.652357,
				0.666035, 0.679238, 0.691978, 0.704266, 0.716114, 0.727533,
				0.738536, 0.749134, 0.759338, 0.769161, 0.778615, 0.787710,
				0.796458, 0.804870, 0.812958, 0.820732, 0.828203, 0.835381,
				0.842276, 0.848898, 0.855256, 0.861361, 0.867222, 0.872846,
				0.878243, 0.883422, 0.888390, 0.893155, 0.897725, 0.902107,
				0.906309, 0.910337, 0.914199, 0.917900, 0.921447, 0.924846,
				0.928102, 0.931222, 0.934211, 0.937074, 0.939816, 0.942442,
				0.944956, 0.947364, 0.949669, 0.951876, 0.953988, 0.956010,
				0.957946, 0.959798, 0.961571, 0.963267, 0.964890, 0.966443,
				0.967928, 0.969350, 0.970709, 0.972010, 0.973253, 0.974443,
				0.975581, 0.976669, 0.977709, 0.978704, 0.979655, 0.980564,
				0.981434};
		for (int n=0; n<octaveResults3degFree.length; n++) {
			assertEquals(octaveResults3degFree[n],
					MathsUtils.chiSquareCdf(n*0.1, 3), 0.000001);
		}

		// Last for 10 degrees of freedom (even, > 2)
		double[] octaveResults10degFree = {0.000000, 0.000000, 0.000000, 
				0.000001, 0.000002, 0.000007, 0.000016, 0.000033, 0.000061,
				0.000106, 0.000172, 0.000266, 0.000394, 0.000565, 0.000786,
				0.001065, 0.001411, 0.001835, 0.002344, 0.002949, 0.003660,
				0.004485, 0.005435, 0.006519, 0.007746, 0.009124, 0.010663,
				0.012370, 0.014253, 0.016320, 0.018576, 0.021028, 0.023682,
				0.026543, 0.029615, 0.032902, 0.036407, 0.040133, 0.044081,
				0.048255, 0.052653, 0.057277, 0.062126, 0.067200, 0.072496,
				0.078014, 0.083751, 0.089703, 0.095869, 0.102243, 0.108822,
				0.115601, 0.122577, 0.129742, 0.137092, 0.144621, 0.152324, 
				0.160193, 0.168223, 0.176406, 0.184737, 0.193207, 0.201811, 
				0.210540, 0.219387, 0.228347, 0.237410, 0.246569, 0.255818, 
				0.265149, 0.274555, 0.284028, 0.293562, 0.303148, 0.312781, 
				0.322452, 0.332156, 0.341886, 0.351635, 0.361396, 0.371163, 
				0.380930, 0.390692, 0.400441, 0.410173, 0.419882, 0.429562, 
				0.439208, 0.448816, 0.458380, 0.467896, 0.477360, 0.486766, 
				0.496111, 0.505391, 0.514602, 0.523741, 0.532804, 0.541788, 
				0.550690, 0.559507};
		for (int n=0; n<octaveResults10degFree.length; n++) {
			assertEquals(octaveResults10degFree[n],
					MathsUtils.chiSquareCdf(n*0.1, 10), 0.000001);
		}
	}
	
	public void testNormalPdf() throws Exception {
		// Check values for x = -4:0.1:4 generated by octave
		double[] expectedPdfMu0Std1 = {0.00013, 0.00020, 0.00029, 0.00042,
				0.00061, 0.00087, 0.00123, 0.00172, 0.00238, 0.00327, 0.00443,
				0.00595, 0.00792, 0.01042, 0.01358, 0.01753, 0.02239, 0.02833,
				0.03547, 0.04398, 0.05399, 0.06562, 0.07895, 0.09405, 0.11092,
				0.12952, 0.14973, 0.17137, 0.19419, 0.21785, 0.24197, 0.26609,
				0.28969, 0.31225, 0.33322, 0.35207, 0.36827, 0.38139, 0.39104,
				0.39695, 0.39894, 0.39695, 0.39104, 0.38139, 0.36827, 0.35207,
				0.33322, 0.31225, 0.28969, 0.26609, 0.24197, 0.21785, 0.19419,
				0.17137, 0.14973, 0.12952, 0.11092, 0.09405, 0.07895, 0.06562,
				0.05399, 0.04398, 0.03547, 0.02833, 0.02239, 0.01753, 0.01358,
				0.01042, 0.00792, 0.00595, 0.00443, 0.00327, 0.00238, 0.00172,
				0.00123, 0.00087, 0.00061, 0.00042, 0.00029, 0.00020, 0.00013};

		// Mean 5.5, std 2.3
		double[] expectedPdfMu5_5Std2_3 = {0.00003, 0.00004, 0.00005, 0.00006,
				0.00007, 0.00008, 0.00010, 0.00011, 0.00014, 0.00016, 0.00019,
				0.00022, 0.00026, 0.00030, 0.00035, 0.00041, 0.00048, 0.00055,
				0.00064, 0.00074, 0.00085, 0.00098, 0.00113, 0.00129, 0.00148,
				0.00169, 0.00193, 0.00219, 0.00249, 0.00283, 0.00320, 0.00361,
				0.00407, 0.00458, 0.00515, 0.00577, 0.00646, 0.00722, 0.00804,
				0.00895, 0.00994, 0.01102, 0.01219, 0.01347, 0.01484, 0.01633,
				0.01793, 0.01965, 0.02150, 0.02347, 0.02558, 0.02783, 0.03021,
				0.03274, 0.03541, 0.03823, 0.04119, 0.04430, 0.04756, 0.05096,
				0.05449, 0.05816, 0.06197, 0.06589, 0.06994, 0.07409, 0.07834,
				0.08267, 0.08708, 0.09156, 0.09608, 0.10063, 0.10520, 0.10978,
				0.11433, 0.11885, 0.12331, 0.12770, 0.13199, 0.13618, 0.14022};

		double x = -4.0;
		for (int xIndex = 0; xIndex < 81; xIndex++) {
			// System.out.printf("%.1f: %.5f (expected %.5f)\n",
			// 		x, MathsUtils.normalPdf(x, mean, stddev), expectedPdfMu0Std1[xIndex]);
			assertEquals(expectedPdfMu0Std1[xIndex], MathsUtils.normalPdf(x, 0.0, 1.0), OCTAVE_RESOLUTION);
			assertEquals(expectedPdfMu5_5Std2_3[xIndex], MathsUtils.normalPdf(x, 5.5, 2.3), OCTAVE_RESOLUTION);
			x += 0.1;
		}
		
	}
	
	public void testNormalCdf() throws Exception {
		// Check values for x = -4:0.1:4 generated by octave
		double[] expectedCdfMu0Std1 = {0.00003, 0.00005, 0.00007, 0.00011, 0.00016,
				0.00023, 0.00034, 0.00048, 0.00069, 0.00097, 0.00135, 0.00187,
				0.00256, 0.00347, 0.00466, 0.00621, 0.00820, 0.01072, 0.01390,
				0.01786, 0.02275, 0.02872, 0.03593, 0.04457, 0.05480, 0.06681,
				0.08076, 0.09680, 0.11507, 0.13567, 0.15866, 0.18406, 0.21186,
				0.24196, 0.27425, 0.30854, 0.34458, 0.38209, 0.42074, 0.46017,
				0.50000, 0.53983, 0.57926, 0.61791, 0.65542, 0.69146, 0.72575,
				0.75804, 0.78814, 0.81594, 0.84134, 0.86433, 0.88493, 0.90320,
				0.91924, 0.93319, 0.94520, 0.95543, 0.96407, 0.97128, 0.97725,
				0.98214, 0.98610, 0.98928, 0.99180, 0.99379, 0.99534, 0.99653,
				0.99744, 0.99813, 0.99865, 0.99903, 0.99931, 0.99952, 0.99966,
				0.99977, 0.99984, 0.99989, 0.99993, 0.99995, 0.99997};

		// Mean 5.5, std 2.3
		double[] expectedCdfMu5_5Std2_3 = {0.00002, 0.00002, 0.00003, 0.00003,
				0.00004, 0.00005, 0.00005, 0.00007, 0.00008, 0.00009, 0.00011,
				0.00013, 0.00015, 0.00018, 0.00021, 0.00025, 0.00030, 0.00035,
				0.00041, 0.00048, 0.00056, 0.00065, 0.00075, 0.00087, 0.00101,
				0.00117, 0.00135, 0.00156, 0.00179, 0.00206, 0.00236, 0.00270,
				0.00308, 0.00351, 0.00400, 0.00454, 0.00516, 0.00584, 0.00660,
				0.00745, 0.00839, 0.00944, 0.01060, 0.01188, 0.01330, 0.01486,
				0.01657, 0.01845, 0.02050, 0.02275, 0.02520, 0.02787, 0.03077,
				0.03392, 0.03733, 0.04101, 0.04498, 0.04925, 0.05384, 0.05877,
				0.06404, 0.06967, 0.07567, 0.08207, 0.08886, 0.09606, 0.10368,
				0.11173, 0.12021, 0.12915, 0.13853, 0.14836, 0.15866, 0.16940,
				0.18061, 0.19227, 0.20438, 0.21693, 0.22991, 0.24332, 0.25714};

		double x = -4.0;
		for (int xIndex = 0; xIndex < 81; xIndex++) {
			// System.out.printf("%.1f: %.5f (expected %.5f)\n",
			// 		x, MathsUtils.normalPdf(x, mean, stddev), expectedPdfMu0Std1[xIndex]);
			assertEquals(expectedCdfMu0Std1[xIndex], MathsUtils.normalCdf(x, 0.0, 1.0), OCTAVE_RESOLUTION);
			assertEquals(expectedCdfMu5_5Std2_3[xIndex], MathsUtils.normalCdf(x, 5.5, 2.3), OCTAVE_RESOLUTION);
			x += 0.1;
		}
		
	}
	
	public void testMultivariateNormalPdf() throws Exception {
		// Testing against data generated using public domain octave code
		//  by Paul Kienzle available at 
		// http://octave.1599824.n4.nabble.com/Multivariate-pdf-of-a-normal-distribution-td1601886.html
		// http://octave.1599824.n4.nabble.com/attachment/1601887/0/mvnpdf.m
		
		double[] means = {1, -1};
		double[][] covariance = {{.9, .4}, {.4, .3}};
		// Evaluated using:
		// for x = -3:0.5:5
		// 	for y = -5:0.5:3
		//	 printf("%.5f, ", mvnpdf([x, y], mu, sigma));
		//  end
		// end
		double[] expectedMvnPdf1 = {0.00000, 0.00000, 0.00000, 0.00001, 0.00005,
				0.00005, 0.00001, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00001, 0.00024, 0.00052, 0.00015, 0.00001, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00001, 0.00052, 0.00289,
				0.00205, 0.00019, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00001, 0.00059, 0.00803, 0.01417, 0.00323, 0.00010, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00033, 0.01129, 0.04944,
				0.02801, 0.00205, 0.00002, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00010, 0.00803, 0.08727, 0.12272, 0.02232, 0.00052, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00001, 0.00289, 0.07789, 0.27187,
				0.12272, 0.00716, 0.00005, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00052, 0.03516, 0.30459, 0.34125, 0.04944, 0.00093, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00005, 0.00803, 0.17257, 0.47987,
				0.17257, 0.00803, 0.00005, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00093, 0.04944, 0.34125, 0.30459, 0.03516, 0.00052, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00005, 0.00716, 0.12272, 0.27187,
				0.07789, 0.00289, 0.00001, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00052, 0.02232, 0.12272, 0.08727, 0.00803, 0.00010, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00002, 0.00205, 0.02801, 0.04944,
				0.01129, 0.00033, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00010, 0.00323, 0.01417, 0.00803, 0.00059, 0.00001, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00019, 0.00205, 0.00289,
				0.00052, 0.00001, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00001, 0.00015, 0.00052, 0.00024, 0.00001, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00001, 0.00005, 0.00005,
				0.00001, 0.00000, 0.00000, 0.00000};

		int expIndex = 0;
		double[] observations = new double[2];
		for (double x1 = -3.0; x1 <= 5.0; x1 += 0.5) {
			for (double x2 = -5.0; x2 <= 3.0; x2 += 0.5) {
				observations[0] = x1;
				observations[1] = x2;
				assertEquals(expectedMvnPdf1[expIndex],
						MathsUtils.normalPdf(observations, means, covariance), OCTAVE_RESOLUTION);
				expIndex++;
			}			
		}

		double[] means2 = {4, 5, -3};
		double[][] covariance2 = {{.9, .4, .25}, {.4, .3, .2}, {.25, .2, .6}};
		// Evaluated using:
		// for x = 2:1:6
		// 	for y = 3:1:7
		//   for z = -5:1:-1
		//	  printf("%.5f, ", mvnpdf([x, y, z], means2, covariance2));
		//   end
		//  end
		// end
		double[] expectedMvnPdf2 = {0.00013, 0.00017, 0.00003, 0.00000, 0.00000,
				0.00393, 0.02507, 0.01871, 0.00163, 0.00002, 0.00001, 0.00033,
				0.00119, 0.00049, 0.00002, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00001,
				0.00001, 0.00000, 0.00000, 0.00000, 0.00705, 0.04084, 0.02764,
				0.00219, 0.00002, 0.00080, 0.02220, 0.07156, 0.02698, 0.00119,
				0.00000, 0.00000, 0.00002, 0.00003, 0.00001, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00082, 0.00433, 0.00266, 0.00019, 0.00000, 0.00383,
				0.09590, 0.28047, 0.09590, 0.00383, 0.00000, 0.00019, 0.00266,
				0.00433, 0.00082, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00001, 0.00003,
				0.00002, 0.00000, 0.00000, 0.00119, 0.02698, 0.07156, 0.02220,
				0.00080, 0.00002, 0.00219, 0.02764, 0.04084, 0.00705, 0.00000,
				0.00000, 0.00000, 0.00001, 0.00001, 0.00000, 0.00000, 0.00000,
				0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
				0.00002, 0.00049, 0.00119, 0.00033, 0.00001, 0.00002, 0.00163,
				0.01871, 0.02507, 0.00393, 0.00000, 0.00000, 0.00003, 0.00017,
				0.00013};

		expIndex = 0;
		observations = new double[3];
		for (double x1 = 2.0; x1 <= 6.0; x1 += 1.0) {
			observations[0] = x1;
			for (double x2 = 3.0; x2 <= 7.0; x2 += 1.0) {
				observations[1] = x2;
				for (double x3 = -5.0; x3 <= -1.0; x3 += 1.0) {
					observations[2] = x3;
					assertEquals(expectedMvnPdf2[expIndex],
							MathsUtils.normalPdf(observations, means2, covariance2),
								OCTAVE_RESOLUTION);
					expIndex++;
				}
			}			
		}
		
		// And now let's test some really random numbers:
		double[] observations3 = {142.1, -10.25, 0.13, 5.13, 5.935};
		double[] means3 = {142.2, -10.3, 0.1, 5, 6};
		double[][] covariance3 = {{1.2, 0.4, 0.6, 0.1, 0.8},
				{0.4, 2, 0.3, 0.2, 0.4}, {0.6, 0.3, 0.9, 0.5, 0.45},
				{0.1, 0.2, 0.5, 1.1, 0.8}, {0.8, 0.4, 0.45, 0.8, 1.3}};
		assertEquals(0.028085,
				MathsUtils.normalPdf(observations3, means3, covariance3),
				0.000001);
	}
	
	public void testDigamma() throws Exception {
		assertEquals(-0.577216, MathsUtils.digamma(1), 0.000001);
		assertEquals(0.42278, MathsUtils.digamma(2), 0.00001);
		assertEquals(2.2518, MathsUtils.digamma(10), 0.0001);
		assertEquals(4.6002, MathsUtils.digamma(100), 0.0001);
		assertEquals(6.9073, MathsUtils.digamma(1000), 0.0001);
		// Test calling for values above the range that we cache
		assertEquals(9.2103, MathsUtils.digamma(10000), 0.0001);
		assertEquals(9.2104, MathsUtils.digamma(10001), 0.0001);
		// Test retrieving a cached value
		assertEquals(5.2958, MathsUtils.digamma(200), 0.0001);

		// And test manually calculated digammas:
		assertEquals(-0.577216, MathsUtils.digammaByDefinition(1), 0.000001);
		assertEquals(0.42278, MathsUtils.digammaByDefinition(2), 0.00001);
		assertEquals(2.2518, MathsUtils.digammaByDefinition(10), 0.0001);
		assertEquals(4.6002, MathsUtils.digammaByDefinition(100), 0.0001);
		assertEquals(6.9073, MathsUtils.digammaByDefinition(1000), 0.0001);
		assertEquals(9.2103, MathsUtils.digammaByDefinition(10000), 0.0001);
		assertEquals(9.2104, MathsUtils.digammaByDefinition(10001), 0.0001);
		assertEquals(5.2958, MathsUtils.digammaByDefinition(200), 0.0001);
	}
}
